home *** CD-ROM | disk | FTP | other *** search
/ Meeting Pearls 1 / Meeting Pearls Vol 1 (1994).iso / installed_progs / gfx / lise2.1 / lise / src / linfit2.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-03-31  |  12.5 KB  |  474 lines
  1. #include <stdio.h>
  2. #ifdef AMIGA
  3. #include <gfxamiga.h>
  4. #else
  5. #include <gfx.h>
  6. #endif
  7. #include <math.h>
  8.  
  9. #define MAXC 10
  10. #define CEXP -1
  11. #define NEXP -2
  12.  
  13. float *x, *y, *y_calc, *resid;
  14. float coef[MAXC], sig[MAXC];
  15. int   nrow, ncol;
  16. float correl_coef;
  17. float *spc, *err, *tim;
  18.  
  19. float frline(stream)
  20. FILE *stream;
  21. {
  22. int i,n;
  23. char c,s[80];
  24.  
  25.    fgets(s,80,stream);
  26.    return(atosf(s));
  27. }
  28.  
  29. square(x, y, a, g, nrow, ncol)
  30.    float  *x[MAXC];
  31.    float  *y;
  32.    float  a[MAXC][MAXC];
  33.    float  g[MAXC];
  34.    int   nrow, ncol;
  35. {
  36. int  i,l,k;
  37.  
  38.    for (k = 0; k <= ncol; k++) {
  39.       for (l = 0; l <= k; l++) {
  40.          a[k][l] = 0.0;
  41.          for (i = 0 ; i <= nrow ; i++) {
  42.             a[k][l] = a[k][l] + x[i][l] * x[i][k];
  43.             if (k != l)   a[l][k] = a[k][l];
  44.          }
  45.       }
  46.       g[k] = 0.0;
  47.       for (i = 0; i <= nrow ; i++)  g[k] = g[k] + y[i] * x[i][k];
  48.    }
  49. }
  50.  
  51. swap(a,b)
  52.    float   *a, *b;
  53. {
  54.    float   hold;
  55.  
  56.    hold = *a; *a = *b;
  57.    *b = hold;
  58. }
  59.  
  60. gaussj(b, y, coef, ncol, error)
  61.    float  b[MAXC][MAXC];
  62.    float  y[MAXC];
  63.    float  coef[MAXC];
  64.    int   ncol;
  65.    int   *error;
  66. {
  67. float w[MAXC][MAXC];
  68. int   index[MAXC][MAXC];
  69. int   irow, icol, n;
  70. int   l1, l, k, j, i;
  71. float determ, pivot, hold, sum, t, ab, big;
  72.  
  73.    *error = FALSE;
  74.    n = ncol;
  75.    for (i = 0; i <= n ; i++) {
  76.       w[i][0] = y[i];
  77.       index[i][2] = 0;
  78.    }
  79.    determ = 1.0;
  80.    for (i = 0; i <= n ; i++) {
  81.       big = 0.0;
  82.       for (j = 0; j <= n ; j++) {
  83.          if (index[j][2] != 1) {
  84.             for (k = 0; k <= n ; k++) {
  85.                if (index[k][2] > 1) {
  86.                   printf("fehler: Matrix singulaer\n");
  87.                   *error = TRUE;
  88.                   return(0);
  89.                }
  90.                if (index[k][2] < 1) {
  91.                   if (((float)fabs((double)b[j][k])) > ((float)big)) {
  92.                      irow = j; icol = k;
  93.                      big = (float)fabs((double)b[j][k]);
  94.                   }
  95.                }
  96.             }
  97.          }
  98.       }
  99.       index[icol][2] = index[icol][2] + 1;
  100.       index[i][0] = irow;
  101.       index[i][1] = icol;
  102.       if (irow != icol) {
  103.          determ = -determ;
  104.          swap(b[irow][0], b[icol][0]);
  105.          swap(w[irow][0], w[icol][0]);
  106.       }
  107.       pivot = b[icol][icol];
  108.       determ = determ * pivot;
  109.       b[icol][icol] = 1.0;
  110.       for (l = 0; l <= n; l++) b[icol][l] = b[icol][l] / pivot;
  111.       w[icol][0] = w[icol][0] / pivot;
  112.       for (l1 = 0; l1 <= n ; l1++) {
  113.          if (l1 != icol) {
  114.             t = b[l1][icol];
  115.             b[l1][icol] = 0.0;
  116.             for (l = 0; l <= n ; l++) b[l1][l] = b[l1][l] - b[icol][l] * t;
  117.             w[l1][0] = w[l1][0] - w[icol][0] * t;
  118.          }
  119.       }
  120.    }
  121.    for (i = 0 ; i <= n ; i++) {
  122.       l = n - i;
  123.       if (index[l][0] != index[l][1]) {
  124.          irow = index[l][0];
  125.          icol = index[l][1];
  126.          for (k = 0; k <= n ; k++) swap(b[k][irow], b[k][icol]);
  127.       }
  128.    }
  129.    for (k = 0; k <= n ; k++) {
  130.       if (index[k][2] != 1) {
  131.          printf("Fehler: Matrix singulaer\n");
  132.          *error = TRUE;
  133.          return(0);
  134.       }
  135.    }
  136.    for (i = 0 ; i <= n ; i++) coef[i] = w[i][0];
  137. }
  138.  
  139. get_data(xarr, yarr)
  140.    float xarr[], yarr[];
  141. {
  142. int   nrow,i,n,x,y,xx,yy,
  143.       leftmargin,
  144.       bottommargin,
  145.       rightmargin,
  146.       topmargin,
  147.       _win_flg;
  148. float tmin,
  149.       tmax,
  150.       ymin,
  151.       ymax,
  152.       real_t,
  153.       real_y,
  154.       yfak,
  155.       tfak,
  156.       xoffset,
  157.       yoffset;
  158. float a,b;
  159. char  s[80];
  160. FILE *fp;
  161.  
  162.    nrow = -1;
  163.    while(1 == 1) {
  164.       printf("%d: X [fp-number | Q<uit> | C<ursor>] x: ",nrow+1); fflush(stdout);
  165.       fgets(s,80,stdin);
  166.       if(toupper(s[0]) == 'Q') break;
  167.       if(toupper(s[0]) == 'C') {
  168.  
  169.          tekopen();
  170.          getxy(&x,&y);
  171.          close(_tek4014);
  172. /* calculating channel and time equivalents */
  173.          fp=fopen(ASHBDRY,"r");
  174.          leftmargin = frline(fp);
  175.          bottommargin = frline(fp);
  176.          rightmargin = frline(fp);
  177.          topmargin = frline(fp);
  178.          _win_flg = frline(fp);
  179.          tmin = frline(fp);
  180.          tmax = frline(fp);
  181.          ymin = frline(fp);
  182.          ymax = frline(fp);
  183.          _tica = frline(fp);
  184.          fclose(fp);
  185.          yfak=(topmargin-bottommargin)/(ymax-ymin);
  186.          tfak=(rightmargin-leftmargin)/(tmax-tmin);
  187.          xoffset= ((- tmin) * tfak) + leftmargin;
  188.          yoffset= ((- ymin) * yfak) + bottommargin;
  189.          a= (x - xoffset)/tfak;
  190.          b= (y - yoffset)/yfak;
  191.          xx= a / _tica;
  192.       }
  193.       if(s[0] <'A') {
  194.          a = atosf(s);
  195.          printf("%d: Y [fp-number]                  y: "); fflush(stdout);
  196.          fgets(s,80,stdin);
  197.          b = atosf(s);
  198.       }
  199.       nrow = nrow +1;
  200.       xarr[nrow] = a;
  201.       yarr[nrow] = b;
  202. printf("x = %f    y = %f\n",a,b);
  203.    }
  204.    return(nrow);
  205. }
  206.  
  207. write_data()
  208. {
  209.     int   i;
  210.  
  211. /* ------------------
  212.    printf("nrow = %d\n\n", nrow);
  213.    printf(" I     X         Y      Y Ber.   Resid\n");
  214.    for (i = 0; i <= nrow ; i++) {
  215.       printf("%d %f %f %f %f\n", i, x[i], y[i], y_calc[i], resid[i]);
  216.    }
  217. ------------------- */
  218.    printf("Koeffizienten  Fehler\n");
  219.    printf("%f     %f    konstanter term\n", coef[0], sig[0]);
  220.    for (i = 1; i <= ncol; i++) {
  221.       printf("%f     %f\n", coef[i], sig[i]);
  222.    }
  223.    printf("korrelationskoeffizient %f\n", correl_coef);
  224. }
  225.  
  226.  
  227. polyfit(x, y, y_calc, resid, coef, sig, nrow, ncol)
  228.    float  *x, *y, *y_calc, *resid, *coef, *sig;
  229.    int   nrow, ncol;
  230. {
  231. float *xmatr[MAXC];
  232. float a[MAXC][MAXC];
  233. float g[MAXC];
  234. int   error;
  235. int   nm;
  236. int   j;
  237. int   i;
  238. float   xi, yi, yc, srs, see, sum_y, sum_y2;
  239.  
  240.    for(i = 0; i < MAXC; i++) {
  241.       xmatr[i] = (float *)calloc(_MAXSPCLEN+2,sizeof(float));
  242.       if(xmatr[i] == NULL) {
  243.          fprintf(stderr,"not enough memory, exiting\n");
  244.          exit(-1);
  245.       }
  246.    }
  247.  
  248.    for (i = 0; i <= nrow; i++) {
  249.       xi = x[i];
  250.       xmatr[i][0] = 1.0;
  251.       for (j = 1; j <= ncol; j++) xmatr[i][j] = xmatr[i][j - 1] * xi;
  252.    }
  253.    square(xmatr, y, a, g, nrow, ncol);
  254.    gaussj(a, g, coef, ncol, &error);
  255.    sum_y = 0.0;
  256.    sum_y2 = 0.0;
  257.    srs = 0.0;
  258.    for (i = 0; i <= nrow; i++) {
  259.       yi = y[i];
  260.       yc = 0.0;
  261.       for (j = 0; j <= ncol; j++) yc = yc + coef[j] * xmatr[i][j];
  262.       y_calc[i] = yc;
  263.       resid[i] = yc - yi;
  264.       srs = srs + (resid[i] * resid[i]);
  265.       sum_y = sum_y + yi;
  266.       sum_y2 = sum_y2 + yi * yi;
  267.    }
  268.    correl_coef = (float)sqrt((double)(1.0 - srs / (sum_y2 - (sum_y * sum_y) / nrow)));
  269.    if (nrow == ncol) {
  270.       nm = 1;
  271.    } else {
  272.       nm = nrow - ncol;
  273.    }
  274.    see = (float) sqrt((double)(srs / (float)nm));
  275.    for (i = 0; i <= ncol; i++) {
  276.        sig[i] = see * (float)sqrt((double)a[i][i]);
  277.    }
  278.    for(i = 0; i < MAXC; i++) free(xmatr[i]);
  279. }
  280.  
  281. help()
  282. {
  283. printf("Linear or Polynomial Fit to Data.\n");
  284. printf("                        Highest possible order: %d\n",MAXC);
  285. printf("LinFIT [-in file] [-o file] [-dif file] [-v] [-poly n]\n");
  286. printf("options:\n");
  287. printf("    -in filename   reads data from spectrum file\n");
  288. printf("    -o  filename   output calculated fitcurve to spectrum file\n");
  289. printf("    -dif filename  diferentiate fitcurve and write data to spectrum\n");
  290. printf("    -v             verbouse: prints out resulting coeffitients\n");
  291. printf("    -poly n        calculate polynom of order n (default is 3)\n");
  292. printf("    -cexp n.m      calculate critical exponent for Tc = n.m\n");
  293. printf("                    Formula: y = a * (Tc - x) ^b\n");
  294. printf("                    if Tc is set to 0, x is assumed to be (Tc - T)\n");
  295. printf("    -nexp n.m      calculate normal exponential curve for Background n.m\n");
  296. printf("                    Formula: y = a * x ^b + n.m\n\n");
  297. }
  298.  
  299. gen_curve(spc,tim,x,y,coef,numchan,ncol,nrow)
  300. float *spc,*tim,*x,*y,*coef;
  301. int numchan,ncol,nrow;
  302. {
  303. float a,b,xx;
  304. int n,m,i;
  305.  
  306. /* --------------- finding maximum and minimum x ------------- */
  307.    a = x[0]; b = a;
  308.    for(i = 0; i <= nrow; i++) {
  309.       if(x[i] < a) a = x[i];
  310.       if(x[i] > b) b = x[i];
  311.    }
  312. /* generating curve */
  313.    for(i = 0 ; i <= numchan ; i++) {
  314.       xx = 1.0; spc[i] = 0.0;
  315.       tim[i] = a + (b - a) * ((float)i)/((float)numchan);
  316.       for(n = 0; n <= ncol; n++) {
  317.          spc[i] = spc[i] + (coef[n] * xx);
  318.          xx = xx * tim[i];
  319.       }
  320.    }
  321. }
  322.  
  323. fit_cexp(x,y,spc,tim,coef,nrow,ncol,tc,numchan)
  324. float *x, *y, *spc, *tim, *coef;
  325. int nrow, ncol;
  326. float tc;
  327. int numchan;
  328. {
  329. int i,n;
  330. float x1,x2,a,b;
  331.    for(i = 0; i <= nrow; i++) {
  332.       x[i] = tim[i]; y[i] = spc[i];
  333.       if(tc != 0.0) x[i] = tc - tim[i];
  334.       if(x[i] <= 0.0) x[i] = 0.00001;
  335.       if(y[i] <= 0.0) y[i] = 0.00001;
  336.       x[i] = (float) log((double) x[i]);
  337.       y[i] = (float) log((double) y[i]);
  338.    }
  339.    polyfit(x, y, y_calc, resid, coef, sig, nrow, 1);
  340.    a = coef[0] ; b = coef[1] ;
  341.    if(b == 0) {
  342.       fprintf(stderr,"sorry, b=0 . exiting.\n");
  343.       exit(0);
  344.    }
  345.    a = (float) exp((double) (a / b));
  346.    x1 = tim[0]; x2 = x1;
  347.    for(i = 0; i <= nrow; i++) {
  348.       if(tim[i] < x1) x1 = tim[i];
  349.       if(tim[i] > x2) x2 = tim[i];
  350.    }
  351. /* generating curve */
  352.    for(i = 0 ; i <= numchan ; i++) {
  353.       tim[i] = x1 + (x2 - x1) * ((float)i)/((float)numchan);
  354.       if(tc != 0.0) {
  355.          spc[i] = a * pow((tc - tim[i]), b);
  356.       } else {
  357.          spc[i] = a * pow(tim[i],b);
  358.       }
  359.    }
  360. }
  361.  
  362. main(argc,argv)
  363. int argc;
  364. char *argv[];
  365. {
  366. int  i,n,m,numchan;
  367. char s[80],z[80];
  368. float xx,fparam,a,b;
  369.  
  370.    checkopt(argc,argv,"-dummy",z);
  371.  
  372.    x      = (float *)calloc(_MAXSPCLEN+2,sizeof(float));
  373.    y      = (float *)calloc(_MAXSPCLEN+2,sizeof(float));
  374.    y_calc = (float *)calloc(_MAXSPCLEN+2,sizeof(float));
  375.    resid  = (float *)calloc(_MAXSPCLEN+2,sizeof(float));
  376.    if(resid == NULL) {
  377.       fprintf(stderr,"not enough memory, exiting...\n");
  378.       exit(-1);
  379.    }
  380.  
  381.    spc    = (float *)calloc(_MAXSPCLEN+2,sizeof(float));
  382.    err    = (float *)calloc(_MAXSPCLEN+2,sizeof(float));
  383.    tim    = (float *)calloc(_MAXSPCLEN+2,sizeof(float));
  384.    if(tim == NULL) {
  385.       fprintf(stderr,"not enough memory, exiting...\n");
  386.       exit(-1);
  387.    }
  388.  
  389.    ncol = 3; numchan = 2048;
  390. /*  ------------------ read order of polynom ------------------ */
  391.    if(checkopt(argc,argv,"-poly",s)) {
  392.       ncol = atoi(s);
  393.       if(ncol > MAXC) {
  394.          printf("you should not try to use more than %d terms\n",MAXC);
  395.          exit(0);
  396.       }
  397.    }
  398.  
  399.    if(checkopt(argc,argv,"-cexp",s)) {
  400.       ncol = CEXP;
  401.       fparam = atosf(s);
  402.    }
  403.  
  404.    if(checkopt(argc,argv,"-nexp",s)) {
  405.       ncol = NEXP;
  406.       fparam = atosf(s);
  407.    }
  408.  
  409. /*  ------------------ read data points ------------------ */
  410.    if(checkopt(argc,argv,"-in",s)) {
  411.       nrow = readspec(s,spc,err,tim,z);
  412.       nrow = nrow - 1;
  413.       for(i = 0; i <= nrow; i++) {
  414.          x[i] = tim[i]; y[i] = spc[i];
  415.       }
  416.    } else {
  417.       nrow = get_data(tim,spc);
  418.       for(i = 0; i <= nrow; i++) {
  419.          x[i] = tim[i]; y[i] = spc[i];
  420.       }
  421.       strcpy(_spc_onam,"");    /* signal "no output file names */
  422.       writespec("userdata",spc,err,nrow+1,'2',"user data");
  423.       writespec("userdata.tim",tim,err,nrow+1,'2',"user data");
  424.       printf("Your data were stored in userdata.spc and userdata.tim\n");
  425.    }
  426.  
  427. /* ------------------ going to start fit ...  ------------------ */
  428. /* --- 1. simple polynomial fit --- */
  429.    if(ncol > 0) {
  430.       polyfit(x, y, y_calc, resid, coef, sig, nrow, ncol);
  431.       gen_curve(spc,tim,x,y,coef,numchan,ncol,nrow);
  432.       if(checkopt(argc,argv,"-v",s)) write_data();
  433.    }
  434. /* --- 2. critical exponent fit --- */
  435.    if(ncol == CEXP) {
  436.       ncol = 1;
  437.       fit_cexp(x,y,spc,tim,coef,nrow,ncol,fparam,numchan);
  438.       a = coef[0] ; b = coef[1] ;
  439.       a = (float) exp((double) (a / b));
  440.       if(checkopt(argc,argv,"-v",s)) printf("Beta = %f  Alpha = %f\n",b,a);
  441.    }
  442. /* --- 3. exponential fit --- */
  443.    if(ncol == NEXP) {
  444. printf("NOT IMPLEMENTED YET !!!!!\n");
  445.       for(i = 0; i <= nrow; i++) {
  446.          x[i] = (float) log((double)(tim[i] - fparam));
  447.          y[i] = (float) log((double) y[i]);
  448.       }
  449.       ncol = 1;
  450.       polyfit(x, y, y_calc, resid, coef, sig, nrow, ncol);
  451.       gen_curve(spc,tim,x,y,coef,numchan,ncol,nrow);
  452.       for(i = 0 ; i <= numchan ; i++) spc[i] = (float) exp((double)spc[i]);
  453.    }
  454.  
  455. /* ------------------ save curve ------------------ */
  456.    if(checkopt(argc,argv,"-o",s)) {
  457.       writespec(s,spc,err,numchan,2,"poly fit curve");
  458.       strcat(s,".tim");
  459.       writespec(s,tim,err,numchan,2,"poly fit curve");
  460.    }
  461. /* ------------------ generate differentiated curve ------------------ */
  462.    if(checkopt(argc,argv,"-dif",s)) {
  463.       for(i = 0 ; i <= (numchan - 1) ; i++) {
  464.          spc[i] = (spc[i] - spc[i+1]) / (tim[i] - tim[i+1]);
  465.       }
  466.       spc[numchan] = spc[numchan-1];
  467.       writespec(s,spc,err,numchan,2,"poly fit differentiated curve");
  468.       strcat(s,".tim");
  469.       writespec(s,tim,err,numchan,2,"poly fit differentiated curve");
  470.    }
  471.    free(x); free(y); free(y_calc); free(resid); free(spc); free(err); free(tim);
  472.    exit(0);
  473. }
  474.